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ABSTRACT 


The  problem  of  estimating  seismicity  and  the  performance  of  a 
system  which  detects  earthquakes  is  formulated  in  such  a  way  that 
maximum  likelihood  estimation  can  be  applied.  The  mean  number  of 
earthquakes  which  occur  in  a  fixed  time  interval  is  assumed  to  be  of  the 
form  exp  |a  -  bmj  where  m  is  magnitude  and  a  and  b  are  constants.  The 
probability  of  detection  as  a  function  of  m  is  taken  to  be  an  error  function 
with  mean  and  variance  \i  and  a.  Procedures  to  obtain  maximum  likelihood 
estimates  of  a,  b,  \l  and  a  are  derived,  discussed,  and  applied  to  experi¬ 
mental  data  to  check  the  relevance  of  the  theoretical  development. 
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I. 


INTRODUCTION 


Suppose  that  a  seismic  station  or  a  network  of  stations  is  operated 
under  a  fixed  set  of  rules  for  the  detection  of  earthquakes  for  a  period  of 
time  of  duration  T.  Suppose  also  that  only  events  which  prove  to  be  earth¬ 
quakes  from  a  preassigned  region  are  retained  and  their  body-wave  magni¬ 
tudes  determined.  The  data  from  such  an  experiment  is  an  integer,  K, 
the  total  number  of  earthquakes  recorded,  and  a  list  of  K  magnitudes, 
m  ,  .  .  .  ,  m  ,  usually  presented  in  the  form  of  a  cumulative  histogram, 
showing  the  log  of  the  number  of  recorded  events  having  magnitude  at  least 
m,  versus  m.  Using  the  available  data  it  should  be  possible  to  determine 
a)  the  seismicity  of  the  region  in  question,  or  b)  the  detection  capability 
of  the  station  or  network  used,  or  c)  both.  In  fact,  one  must  almost  always 
determine  both  seismicity  and  detection  performance,  although  the  chief 
objective  may  be  to  find  only  one  of  these  quantities,  since  the  other  is 
rarely  well  enough  known  in  advance  to  be  modeled  accurately. 

The  aim  of  this  note  is  to  obtain  simple  rules  for  the  estimation  of 
seismicity  and  system  performance  from  the  data  of  such  an  experiment, 
and  to  get  some  idea  of  the  statistical  validity  of  these  estimates.  In  order 
to  reduce  this  problem  to  the  estimation  of  specific  parameters  from  the 
data  (so  that  standard  statistical  estimation  techniques  may  be  applied),  it 
is  necessary  to  make  specific  assumptions  about  the  random  character  of 
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seismicity  and  the  detection  of  weak  events  in  background  noise  by  the 
seismic  network.  In  order  to  push  the  mathematics  through  we  make  use 
of  the  well-known  and  tractable  laws  of  Poisson  and  Gauss,  and  our  results 
are  only  as  valid  as  our  underlying  assumptions,  although  the  method  is 
quite  general. 

Specifically,  we  assume  that  earthquakes  occur  as  a  Poisson  process, 

with  rate  depending  upon  magnitude  and  other  parameters.  Since  in  the 

type  of  experiment  of  interest  here  one  is  concerned  only  with  total  numbers 

of  events  which  occur  over  a  fixed  period  of  time,  we  need  make  no  claim 

that  the  Poisson  process  accurately  represents  the  occurrence  of  events  as 

a  stochastic  process,  but  only  that  the  numbers  of  events  recorded  having 

magnitudes  (and  other  parameters)  in  given  ranges  are  Poisson  variables. 

Let  N  be  the  mean  number  of  events  which  occur  during  T  having  magni- 
m 

tude  at  least  m.  In  addition  to  the  Poisson  assumption  further  we  assume 
a  linear  relationship, 

log  N  =  a  -  bm  ,  ( 1 ) 

m 

between  log  N  (natural  logarithm)  and  magnitude.  If  base-ten  logarithms 
are  used,  a  and  b  become  a*  =  a/log  10  and  b*  =  b/log  10.  Finally,  the 
probability  n(m),  that  the  seismic  system  will  detect  an  event  having 
magnitude  m  (which  will  be  more  precisely  defined  in  Section  HI)  is  assumed 
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to  be  an  error  function: 


2 

n(m)  =  {Zn  a  )  I  exp  -  - - -  dm'  (2) 

2a 

In  terms  of  these  models,  the  general  problem  is  the  estimation  of 

the  four  parameters  a,  b,  \Jl,  and  a  from  the  K  +  1  data  values 

K,  m  ,  .  .  .  ,  m  .  The  quantities  a  and  b  define  the  seismicity  of  the  region, 
1  K 

and  p  and  a  the  detection  capabilities  of  the  system.  In  Section II,  under 
the  simplifying  assumption  that  a  =  0  and  that  \Jl  is  known,  we  discuss  some 
of  the  implications  of  the  Poisson  description  of  natural  seismicity.  In 
Section  3  we  formulate  and  solve  the  estimation  problem  in  the  general 
case,  and  in  Section IVwe  illustrate  the  technique  by  some  application  to 
real  data. 
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IL  ESTIMATING  SEISMICITY 


We  assume  that  an  earthquake  may  be  described  by  its  magnitude,  m, 
and  a  finite  set  of  other  parameters,  symbolized  by  the  vector  a,  with 
sufficient  accuracy  so  that  the  a  priori  probability  of  its  detection  by  a 
given  station  or  network  may  be  determined  if  m  and  a  are  known.  The 
network  must  be  specified  in  considerable  detail,  of  course,  for  such  a 
probability  function  to  be  meaningful.  This  specification  will  include  the 
nature  of  the  sensors  and  recording  equipment  used,  station  noise  levels 
(and  probability  distributions)  as  well  as  the  multi  -  station  data  processing 
algorithm  which  is  used  for  event  detection.  The  earthquake  parameters, 
will  include  latitude  and  longitude  of  epicenter,  depth,  fault  plane 
orientation  parameters,  etc.  At  the  very  least,  knowledge  of  the  network 
and  specification  of  the  parameters  m  and  a  should  permit  the  computation 
of  peak  signal-to-noise  ratios  at  each  of  the  stations  of  the  system. 

Suppose  that  A  stands  for  a  set  of  possible  values  of  the  (vector) 
parameter  a,  and  that  A^  is  the  set  of  a-values  which  might  exist  for  the 
region  under  study.  The  Poisson  '*rate  function’*  for  events  of  magnitude  m 
and  parameter  value  a  will  be  denoted  by  V(m,Q'),  and  we  assume  that  the 
probability  of  occurrence  of  an  earthquake,  in  any  infinitesimal  interval  dt, 
having  a  magnitude  in  the  range  m^^  to  m^  and  a- parameter  in  some  set  A  is 


(3) 
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Furthermore,  according  to  the  Poisson  hypothesis,  the  occurrences  which 
take  place  in  non-overlapping  time  intervals,  or  with  different  magnitudes 
or  a-values,  are  independent  events. 

The  rate  of  occurrence  of  events  specified  by  magnitude  alone  will  be 
called  v(m): 


v(m) 


=  v(m,  £_)  d  a_ 

00 


(4) 


Thus,  the  probability  of  occurrence  of  an  event,  having  magnitude  at  least 
m,  in  time  dt  is 


v(m* )  dm*  dt 


and  the  expected  number  of  such  events  in  a  finite  interval  of  duration  T  is 

called  N  : 

m 


N  =  T  f  v(m> )  dm* 
m 


We  assume  that  N  has  the  form 
m 


N  =  e 
m 


a-bm 


whence,  by  differentiation,  we  obtain 

/  X  rr,  1  a-bm 
v(m)  T  =  b  e 


(5) 


(6) 


(7) 


The  number  of  events,  having  magnitudes  in  the  range  ^m 
which  actually  occur  during  T  is  a  random  variable,  say,  Nj2*  The 


5 


(8) 


expected  value  of  is 


N 


/ 

rr.  r  2  ,  ,  ,  a  /  -bmi  - bm-,  \ 

=  T  v(m)  dm  =  e  ye  ^  -  e 


and  the  probability  that  actually  has  the  value  k  is  given  by  the  Poisson 
distribution: 


Prob  =  k} 


(9) 


Equation  (9)  is  the  result  of  our  basic  assumptions,  and  it  also  follows  that 
if  we  keep  track  of  the  numbers  of  events  which  occur  during  T  having 
magnitudes  in  a  whole  set  of  non-overlapping  magnitude  intervals,  then 
these  numbers  are  independent  random  variables,  each  having  a  distribu¬ 
tion  of  the  form  (9)  with  an  appropriate  mean  value.  From  this  fact  it  is 
possible  to  compute  the  conditional  probability  that  value  k, 

given  that  the  total  number  of  events  having  magnitude  at  least  jj  (call  this 
random  variable  N  )  is  K,  where  K  >  k  and  ll  is  smaller  than  m,  .  This 

\ji  —  1 

probability  is 


Prob  |n,^ 

I  12 


(10) 


According  to  (6), 


N 


exp  (a-b|a). 
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From  the  probability  distribution  (10),  which  is  simply  the  familiar 


binomial  distribution,  giving  the  chance  of  k  successes  out  of  K  trials, 
with  a  success  probability  of  N,^/N  ,  we  can  easily  compute  the  conditional 

12  fi 

mean  and  variance  of  ^^2'  given  that  N  =  K.  The  conditional  mean, 

EK(N,^)is 


N 

N  =  KJ-  =  K  =: 


12 


N 


and  the  standard  deviation,  a  (N  )  (given  that  N  =  K),  is 

K  12 


(11) 


(12) 


We  illustrate  some  implications  of  the  formulas  by  discussing  the 
interpretation  of  a  hypothetical  experiment  designed  to  study  seismicity  for 
large  events.  We  assume  that  one  is  interested  only  in  the  rate  of  occur¬ 
rence  of  clearly  detectable  events  and  that  he  is  willing  to  assume  that  all 
events  of  magnitude  at  least  [i  are  recorded  by  the  system.  The  detection- 
system  parameters  are  then  removed  from  the  problem  and  it  remains  to 
estimate  the  seismicity  parameters,  a  and  b,  from  the  data.  The  remainder 
of  this  section  is  devoted  to  this  problem,  and  the  more  general  case  in 
which  the  detection  capability  of  the  network  is  also  modeled  is  treated  in 
later  sections.  As  before,  we  assume  that  K  events  are  recorded,  with 
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magnitude sm.,  i=  all  at  least  equal  to  [i.  Let  N  stand  for  the 


m 


number  of  recorded  events  having  magnitude  at  least  m.  Then  N  ,  the 


m 


unconditioned  mean,  is  N  =  exp  (a-bm),  and  the  conditional  mean,  given 


m 


that  N  =  K,  is  given  by  (11): 


E^(N  )  =  =  K 


(13) 


The  usual  display  of  such  experimental  data  is  a  plot  of  log  N  versus  m. 


m 


In  such  a  plot,  log  (E  (N  ))  will  appear  as  a  straight  line  of  slope  (-b).  If 

K  m 

base-ten  logarithms  are  used,  we  rewrite  (13)  as 


E  (N  )  =  K  10 
K  m 


-b*(m-|a) 


(14) 


where  b*  =  b/log  10  =  0,434Z9...b. 

As  m  increases,  the  expected  number  of  events  decreases  rapidly, 

and  with  it,  the  ratio  of  standard  deviation  (formula  (12))  to  mean  increases. 

1/2 

In  fact,  for  very  small  values  of  E^,(N  ),  a^,(N  )  ?^  (E^,(N  ))  .  Since 

K  m  K  m  K  m 

the  conditional  probability  distribution  for  (given  =  K)  is  the 

binomial  law  (10),  we  can  establish  intervals ,  about  the  value  of  E  (N  ), 

i\  m 

which  contain  the  value  of  N  with  fixed  probability.  This  is  illustrated  in 

m 

Figures  1  and  2  for  K  =  100  and  K  =  1000,  where  the  solid  lines  are  plots  of 


E^,(N  )  versus  m  (assuming  that  b*  =  1.0)  and  the  dashed  lines  contain  the 

Km 
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values  of  N  with  50%  probability  (i.e.  ,  they  mark  the  points  outside  of 
m 

which  the  tails  of  the  distribution  contain  25%  probability)  and  90%  probabi¬ 
lity  (i.  e,  ,  5%  tails),  as  labeled. 

These  illustrations  show  that  it  is  difficult  to  fit  a  straight  line  to 
data,  presented  graphically  in  this  form,  and  put  the  proper  statistical 
weight  to  the  various  portions  of  the  curve.  The  determination  of  a  straight 
line  fit  by  simple  leas t- squa res  has  also  been  criticized  on  this  basis  (2). 
The  problem  is  seriously  aggravated  in  the  usual  case,  where  the  cumula¬ 
tive  histogram  begins  to  level  off,  going  towards  lower  magnitudes,  because 
the  detection  system  begins  to  miss  events,  just  as  the  number  of  recorded 
events  is  beginning  to  give  a  statistically  reliable  picture  of  the  seismicity 
law.  For  these  reasons,  statistical  estimation  theory  is  applied  to  this 
problem,  using  the  maximum  likelihood  method  for  the  determination  of 
unknown  parameters. 

From  the  independence  of  the  numbers  of  events  in  different  magni¬ 
tude  ranges,  it  follows  that  the  conditional  probability,  given  =  K,  that 
the  K  recorded  events  actually  have  magnitudes  in  the  intervals  m.  to 


m.  +  dm. ,  i  =  1,...K,  is 
1  1 


k!  .n, 
1  =  1 


v(m.)  T  dm.  i  K  i 

^  i  i  .,i,K  bjjK  ^  -bm^ 


N 


=  K!b^^  e  ^  .He  dm. 

1=  1  1 


li 


(15) 
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We  assume  that  the  observed  values  of  magnitude,  m  .  .  .m  ,  are  arranged 

1  K 

in  non-decreasing  order  (not  in  the  order  of  occurrence).  Thus 

P  (m  .  .  •  is  a  probability  density  function,  defined  and  normalized 
K  1  K 

over  the  portion  of  K-space  given  by 

-•••  -"k  • 

The  corresponding  likelihood  function,  A  is  given  by  log  P  ,  and  the 

K  K 

maximum  likelihood  estimator  of  b  is  the  value  which  maximizes  A  We 

K . 

denote  by  (m)  the  average  observed  magnitude: 


<m> 


i. 


K  1  =  1 


(16) 


and  we  compute 

Aj^(m^ ,  .  .  .  ,  m^  :  b)  =  K 


log  b  +  (p-  <m))  b 


+  log  K! 


The  value,  b,  which  maximizes  A  is  obviously 

K 


b  =  (<m>  -  p) 


(17) 


Because  of  the  form  of  (17),  it  is  more  convenient  to  parametrize  the 
seismicity  by  an  inverse  parameter, 


m  =  1/b 
e 


(18) 


so  that  the  mean  seismicity  curve  has  the  form 


E  (N  )  =  e 
K  m 


a-(m/m0) 


(19) 


Then  the  estimator  of  m  corresponding  to  (17)  is 

e 
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=  <m>  -  , 


(20) 


A 

m 

e 

Equation  (20)  implies  that  the  seismicity  parameter,  m^,  is  chosen 
so  that  the  expected  number  of  events  of  magnitude  at  least  equal  to  the 
sample  mean  (m)  is  just  (1/e)  times  the  number  actually  contained  in  the 
sample  (which  contains  only  events  of  magnitude  at  least  |i).  If  one  had  a 
really  large  sample  of  events  covering  a  large  range  of  magnitudes,  it 
would  be  possible  to  test  the  assumption  of  exponential  seismicity,  as 
expressed  by  (19),  by  treating  the  data  as  a  series  of  experiments,  each 
containing  the  data  of  its  predecessors,  and  each  having  a  smaller  value  of 

One  can  get  some  idea  of  the  accuracy  of  the  estimator  (20)  by  compu¬ 
ting  its  mean  and  standard  deviation,  conditioned  on  N  =  K.  By  direct 
integration,  using  the  conditional  probability  density  (15)  we  find  that 

E^((m)  )  =  b  e^^  ^  ^  mdm=|jL+l/b 
or 

E  (m  )  =  1/b  =  m  .  (21) 

K  e  e 

In  other  words,  our  estimator  of  m  is  unbiased.  We  also  find 

e 

E  (<m>^)  =  (ij  +  m  )^  +  (m  /K)  , 

K  e  e 

so  that  the  standard  deviation  of  m  (given  N  =  K)  is 
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(22) 


CTT^(m  )  =  m  /K 
K  e  e 


1/2 


If  we  put  b'  =  1/m  ,  so  that  (19)  is  the  same  as 
d 

E  (N  )  = 

K  m 

then  m,  =  (log  10)m  and  m,  =  (log  10)  •  m  ,  Both  estimators  have  the 
d  e  d  e 

same  relative  accuracy: 


OK(me> 


=  K 


-1/2 


(23) 


Thus,  a  1  00-ear thquake  sample  should  provide  10%  accuracy  in  the 
measurement  of  seismicity. 

So  far  we  have  ignored  the  parameter  a,  since  it  drops  out  of  the 
conditional  probability  density  from  which  b  was  estimated.  But  the  fact 
that  K  events  were  seen,  with  magnitude  at  least  |a,  is  itself  a  measure  of 
total  seismicity,  and  hence  can  be  used  to  determine  a.  According  to  our 
basic  model,  the  probability  that  N  is  equal  to  K  is  given  by  (9): 


Prob  (n 
1  U 


(N  )^ 
K! 


(24) 


where  N  =  exp  (a-bja),  and  b  has  already  been  estimated.  The  maximum 
likelihood  estimator  of  a  is  that  value,  a,  which  maximizes  (24),  for  fixed 
K.  But  the  function  X  exp  (-X)  is  a  maximum  at  X  =  K,  hence  the 
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estimate,  a,  must  make  N  =  K,  i.  e.  , 


a  =  b  +  log  K 


(25) 


If  we  approach  the  problem  in  terms  of  the  joint  estimation  of  a  and  b,  we 

are  led  immediately  to  the  same  results,  as  will  be  seen  in  the  next 

section,  when  a  more  general  problem  is  solved. 

It  should  be  observed  that  equation  (15)  lends  itself  to  a  slightly  differ 

ent  interpretation.  Instead  of  thinking  of  a  single  experiment  which  yields 

the  magnitude  values  m  ,  m  ,  .  .  .  ,  m  (listed  in  non-decreasing  order ),  we 

12  K 

may  think  of  the  result  as  the  outcome  of  K  independent  trials  of  a  single 
experiment,  each  of  which  yields  a  single  value  of  magnitude.  The  probabi 
lity  density  for  magnitude,  now  treated  as  a  random  variable,  is 


f(m) 


v(m)T 

SLl  - 

N 

U 


be 


which  is  normalized  in  the  interval  |i  <.  “  •  If  is  the  magnitude  of  the 
th 

i  event,  then  the  probability  density  for  the  sequence  m  ,  m  ,  ...,m  is, 

12  i\ 

of  course, 

iSl  . 

Finally,  the  probability  density  for  obtaining  the  values  m^,  m^,  .  .  .  ,  m^ 

(listed  in  non-decreasing  order)  will  be  the  sum  over  all  possible  orders  at 

occurrence,  i.  e.  ,  K!  F  (m  ,  .  .  .  ,  m  ),  which  is  exactly  F  (m  ,  .  .  .  ,  m  ), 

K  1  K  i\  1  x\ 
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as  given  by  (15).  In  the  literature  (1  )(Z  )(4)(5)(6)  the  problem  has  been 
treated  from  this  second  point  of  view,  starting  with  f(m)  and  F(m  ,  .  .  .  ,m  ) 
as  the  basic  probability  laws.  In  this  way  the  estimator  (17)  has  been 
found,  and  formulas  (Zl)  and  (ZZ)  derived.  In  fact,  the  probability  density 

A, 

of  b  has  been  obtained  and  confidence  intervals  established  for  this  estima¬ 
tor  by  Aki  and  Utsu. 


14 


III.  ESTIMATING  NETWORK  PERFORMANCE 


In  Section  II  we  introduced  the  '*rate  function,'*  which 

describes  the  relative  rate  of  occurrence  of  earthquakes  in  terms  of  magni¬ 
tude  and  the  vector  parameter  a.  The  ratio  of  this  rate  to  the  integrated 
rate  v(m)  will  be  denoted  by 


((a  I  m) 


V(m,  a) 
V(m) 


(26) 


Because  of  the  meaning  of  v(m,a’)  dm  do'  dt  and  v(m)  dm  dt  as  probabilities, 
f(a’|m)  itself  is  the  conditional  probability  density  that,  given  the  occurrence 
of  an  event  with  magnitude  m,  the  event  had  a  parameter -value  a.  We  have 
assumed  that  a  and  m  together  characterize  an  event  well  enough  so  that  its 
probability  of  detection  by  the  network  can  be  specified.  Let  this  detection 
probability  be  P(m,Q').  Then  the  probability  that,  in  time  dt,  an  event  in 
the  magnitude  range  m  to  m  +  dm  will  occur  and  be  detected  is 


J  P(m,Q')  vim,^)  da  dm  dt 
*Aoo 

=  v(m)  dm  dt  J  P(m,  a)  f(a|  m)  da 
Ac»  ” 

We  define  the  integral  which  appears  here  to  be  n(m): 


n{m)  =  J  P(m,  O’)  f (o’ 1  m)  dcf 


(27) 


This  quantity  is  clearly  the  probability  that  the  network  detects  an  event 
of  magnitude  m,  knowing  only  m,  but  taking  account  of  the  different  ways 
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this  can  happen  (different  a- values)  and  weighting  them  according  to  the 
detailed  seismicity  of  the  region,  as  described  by  f(a|m).  In  general,  we 
do  not  know  f(a|m),  and  cannot  readily  estimate  it  from  the  data,  hence  we 
must  work  directly  with  Il(m),  the  effective  detection  probability,  and  we 
will  make  simple  assumptions  about  the  form  of  this  function,  in  spite  of 
its  dependence  on  the  nature  of  the  region  studies.  In  terms  of  II(m),  the 
rate  of  recorded  events  over  an  interval  of  duration  T,  having  magnitudes 
in  the  range  m<  is 

^12  ”  ^  'Ll  =  he  ^  n(m)  e  ^dm  (28) 

The  effective  detection  probability,  n(m),  must  contain  parameters 
characteristic  of  the  system  which  can  be  estimated  from  the  data.  For  one 
of  these  parameters  we  choose  the  magnitude,  |ji,  for  which  the  effective 
detection  probability  is  50%,  i.  e.  ,  we  put 

Il(m)  =  II  (m-|a)  (29) 

o 

where  11  (m)  is  a  function  which  increases  from  zero  (at  m-*  -00)  to  unity 
o 

(m-  +  00),  and  goes  through  the  value  1/2  at  m  =  0.  Later,  we  shall 
specialize  Il^(m)  to  be  an  error  function,  but  it  is  convenient  to  carry  out 
the  analysis  using  form  (29)  as  far  as  the  estimation  of  a,  b,  and  jj.  At  that 
point  we  shall  introduce  the  error-function  hypothesis  for  n^(m)  (the  reader 
may  be  moved  to  try  other  models)  and  estimate  also  the  variance 
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parameter,  a  (see  equation  2). 


Thus,  the  expected  number  of  recorded  events  in  the  magnitude  range 
m^  S  rn  ^  (formula  28)  becomes 


Ni^  =  be 


a-b^  prr>2-M 


^  1 


-bm 


n  (m)  e  dm 
o 


(30) 


From  (30),  we  see  that  the  expected  number  of  recorded  events  of  magni¬ 
tude  at  least  m  is 


N 

m 


b  e 


■b^  p® 
mi  • 


n  (m*: 

o 


“bm* 


dm* 


(31) 


and  the  expected  total  number  of  recorded  events  is 

N  =  b  1(b)  ,  (32) 

where 

1(b)  =  r  n  (m)  e  dm  .  (33) 

*^-00  O 

It  is  assumed  that  11  (m)  vanishes  rapidly  enough,  as  m-^  -oo,  to  guarantee 
the  convergence  of  this  integral. 

We  can  now  compute  the  likelihood  function  for  the  problem  and 
formulate  the  estimation  of  all  the  unknown  parameters.  As  before,  the 
experiment  consists  of  observation  of  a  given  region  for  a  time  of  duration 
T,  using  a  fixed  detection  network,  and  this  time  K  represents  the  total 
number  of  events  recorded,  regardless  of  magnitude,  while  the  actual 
magnitudes  observed  are  m.,  i  =  1,...,K.  For  a  given  set  of  parameters , 
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a,  b,  |a  (and  others  implicit  in  11  (m)),  the  probability  of  recording  just 
K  events  in  time  T  is  the  Poisson  law 


Pr  ob 


{n  =  k} 


(N)^  -N 


(34) 


with  N  given  by  (32).  If  we  define  K  non-overlapping  magnitude  intervals 

““  th 

and  if  N,  is  the  expected  number  of  events  in  the  i  interval  (given  by 

formula  (30)  with  appropriate  values  of  m^  and  ni^)*  then  the  conditional 

probability,  given  K  detections,  of  obtaining  just  one  event  in  each  of  these 

K  magnitude  intervals  is 


K  _  _ 

Kl  .n,  (N./N) 

1=1  1 


(35) 


.  th 


When  the  magnitude  intervals  are  made  infinitesimal,  so  that  the  i  ranges 

from  m  to  m  +  dm  ,  we  have 
111 

N.  =  T  n(m  )  v(m.)  dm. 

1  111 

and  the  total  probability  of  recording  K  events  with  magnitudes  in  these  K 
intervals  is 


P(K,  m  ,  .  .  .  ,  m_)  dm  .  .  .  dm 
1  ]\  1  Jx 


where 


-N  K 

P(K,  m, . =  e  .n  T  n(m. )  v(m. ) 

1  K  1=1  11 


-N  K 
=  e  be 


|a-b<m>j 


K 


iil  n^(ni.-b) 


(36) 
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In  (36)  we  have  again  made  use  of  the  symbol  (m)  for  the  average  of  the 

observed  magnitudes,  as  given  by  16),  and  N  is  given  explicitly  by  (32)* 

P(K,  m^,  .  ,  .  absolute  probability,  not  a  conditional  probability 

like  P  (m  .  .  .  m  ).  The  expression  P(K,  m  .  .  .  m  )  is  also  defined  for 
ix  i  lx  IK 

magnitude  variables  satisfying  m  <  m  <.  .  .  <  m  .  Finally,  the  likelihood 

12  K 

function  is 


A(K,  m  ^ ,  .  ,  ,  ,  ; 


=  -N  +  K 


j^a  -  b(nn) 


a,  b,  la)  =  log  P(K,  . m^) 

+  log  bj  +  log  n^(nn.-ia) 


(37) 


The  maximum -likelihood  estimators  of  a,  b  and  la  are  obtained  by 
maximizing  A  over  these  parameters.  Differentiating  with  respect  to  a, 
we  require  that 


SA  3n  _  ^ 

3a  "  "Sa  +  ^  ° 

But,  according  to  (32), 

SN  - 

3a  ' 

hence,  the  estimator,  a,  is  that  value  of  a  for  which  (with  given  b,|i) 

N  =  K.  This  result  is  already  interesting  in  that  it  states  that  the  overall 
seismicity,  which  if  governed  by  the  parameter  a,  must  be  chosen  so  that 
the  expected  total  number  of  recorded  events  is  equal  to  the  actual  number 
seen.  It  is  easily  seen  that  (34)  is  maximized  when  N  =  K,  hence  estimating 
the  mean  of  a  sample  Poisson  distribution  from  a  single  measurement 
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produces  the  same  result  we  have  found  here  for  the  more  complex  situation 


of  a  Poisson  process  with  parameters. 

From  (32)  we  can  write 

log  N  =  a  -  b|j  +  log  b  +  log  1(b)  , 

hence,  the  estimator  which  equates  N  to  K  is 

a  =  log  K  +  b|a  -  log  b  -  log  1(b)  .  (38) 

Substituting  this  value  into  (37)  we  obtain 

,  .  .  .  ,  m^  ;  a,  b,  |a)  =  K  ^log  K-l+b|ji-b<m>  -  log  1(b)  j 
K 

+  ill  (m.-n)  .  (39) 


The  equations  for  the  estimation  of  b  and  ^  are  obtained  by  setting  the 

partial  derivatives  of  A  ,  given  by  (39),  with  respect  to  these  parameters 

K 

equal  to  zero: 


|i  -  <m>  - 


I'(b) 

Kb) 


=  0 


1 


K 

-  S  iii 


n  '(nn.-i-i) 

O  1 

n  (m.-n) 

o  1 


=  0 


(40a) 

(40b) 


In  equations  (40)  we  have  used  the  prime  to  denote  the  derivative  of  a  func¬ 
tion  with  respect  to  its  argument.  Thus, 

=  .  .  r  n  (41) 

db  -00  o 
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(we  assume  that  11  (m)  vanishes  rapidly  enough,  as  m-»  -oo,  so  that  the  first 
several  derivatives  of  1(b)  exist). 

In  order  to  proceed,  we  must  now  make  an  assumption  for  11  (m), 

o 

perhaps  containing  further  parameters.  We  shall  then  have  explicit  equa¬ 
tions  for  [x  and  b,  and,  by  differentiating  A  with  respect  to  the  new  para¬ 
meters,  we  get  equations  for  them  as  well.  In  passing  we  note  that  the 
case  treated  in  Section  II  can  be  retrieved  by  treating  |a  as  known  (hence 
ignoring  (40b))  and  putting  n^(m)  equal  to  zero  for  negative  m,  unity  for 
positive  m.  We  also  assume  that  each  m.  <  |a,  since  the  contrary  would  be 
contradictory  and  would  also  make  P  vanish.  Then  1(b)  becomes  simply 
1/b,  and  (40a)  immediately  reproduces  (17). 

As  promised  in  the  Introduction,  we  shall  assume  an  error-function 
law  for  n^(m).  We  shall  ultimately  show  by  example  in  Section  IV  that  it 
can  provide  a  reasonable  representation  of  experimental  results.  Define 

erf(x)  =  (Ztt  )  exp  (-y^/2)  dy  (42) 

—00 

and  take 

n  (m,  a)  =  erf{m/a)  =  2tt  exp{-y^ /2a^  )  dy  ,  (43) 

O  -00 

With  this  assumption  we  find  (integrating  by  parts) 

,  1  _  ,  bn  (m.a) 

-bm  1  -bm  o 

1(b)  =  I(b,  a)  =  J_^  n^{m,  a)  e  dm  =  -J_^  e  — ^ - dm  . 
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But 


an  (m, a)  ?  _i /  ?  7 

— -  =  (2tt  a  )  exp  (-m  /2a)  , 


the  Gauss  function,  hence, 


T/U  \  /O  2-1/2 
I(b,  a)  =  (Ztt  b  a  )  J  e 


m 


2a 


+  bm 


1  ^ - 

dm  =  p  e  ^  (44) 

b 


2 

b  a 


It  follows  from  this  that 


k2  2 

?  1  ^ 

I'(b)  =  (o'"  -  — )  e“2— 


Using  the  error  function  for  II  ,  equation  (40a))  for  the  determination 


of  b  is  simply 


t  2  1  /  \ 

ba  =(a-<m) 

D 


(45a) 


Again,  we  replace  b  by  its  inverse,  m  ,  which  parameter  is  estimated  by 


m  ,  determined  from 
e 

2 

m  -  =  (m)  -  |i  .  (45b) 

e  m 

e 


In  the  limit  a 0,  our  function  11  (m,a)  reverts  to  the  simple  zero-or-one 


probability  which  describes  the  case  treated  in  Section  II  and  (45b) 


obviously  reduces  to  equation  (20)  of  that  section.  In  general,  the  solution 


of  (45b),  which  reduces  to  the  proper  limit  as  a  -♦  0,  is 
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=  (m>  -  ^  + 


1  +  \/l  + 


4a 


4 


<m>  -  ij  ,/  \  ,3 

(<m>  -  ij) 


+  . 


(46) 


In  expression  (46)  the  parameters  of  the  detection  system,  \i  and  a, 
are  yet  to  be  estimated.  When  their  estimates  are  found,  they  are  sub¬ 
stituted  in  (46)  and  (38)  to  provide  the  final  estimates  of  the  seismicity 
parameters.  However,  if  one  felt  that  the  detection  system  were  ade¬ 
quately  represented  by  n(m),  as  given  by  (29)  and  (43)  with  known  values 

of  and  a,  then  (38)  and  (46)  would  be  final  estimates  of  a  and  m  .  It  is 

e 

logical  in  that  case  to  ask  about  the  statistical  properties  of  m  as  an 

estimator  of  m  .  However,  the  form  of  (46)  makes  a  direct  evaluation  of 
e 

the  conditional  mean  and  variance  of  m  impossible.  The  best  we  can  do  is 

e 

evaluate  the  mean  and  variance  of  (m),  as  we  did  in  Section  II  in  a  simpler 
problem,  and  relate  these  quantities  indirectly  to  m^,  or  b,  through  (40a). 
These  calculations  are  similar  to  those  leading  to  (21)  and  (22),  and  we 
obtain 


and 


)  = 


b  e^  r  e  n  (m-|a)m  dm 

'J-oo  o 

a  -bm 

b  e  f  e  n  (m-^)  dm 

'J-oo  o 


1(b)  ^■l(b) 
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Thus 


1  JI"(b) 
K  ll(b) 


I'(b)  Z-i 
1(b)  J 


-  ^)  =  - 


I'(b) 

1(b) 


(47a) 


and 


^)  =  K 


-1/2 


ri"(b) 

ll(b) 


(47b) 


Let  us  put 


A 

Then  (40a)  says  that  the  estimator,  b,  is  determined  by 
-$(b)  =  <m>  -  \i 
and  (47a)  says,  in  turn,  that 

Ej^  ($(b))  =  $(b)  ,  (49a) 

i.e.  ,  $(b)  is  an  unbiased  estimator  of  $(b),  if  |a  and  a  are  known.  Moreover, 
equation  (47b)  is  equivalent  to 


($(b))  =  [$'(b)/K]^^^  (49b) 

A 

Equations  (49)  together  show  that,  as  K ,  the  ’'estimator'*  $(b)  converges 
(in  the  mean)  to  $(b),  and  in  this  sense  b  must  become  a  better  and  better 
estimate  of  b. 

To  continue  with  the  general  problem,  we  now  have  explicit  formulas 
for  a  3-nd  b  in  terms  of  [j  and  a.  In  order  to  determine  |i  and  a  we  must 
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minimize  A(K, 


m  ,...,m  ;a,  b,  \Jl,  g)  over  these  parameters.  Since 

1  K 

equation  (40b),  equivalent  to  BA/S|a  =  0,  and  a  similar  equation,  SA/5a  =  0, 
cannot  be  solved  explicitly,  one  has  to  proceed  numerically,  and  for  this 
purpose  it  is  simpler  to  work  with  the  likelihood  function  itself,  seeking  a 
minimum  by  some  search  procedure.  According  to  (45a), 

b(u  -  <m>)  =  (ba)^  -  1 
and  also,  from  (44) 

log  I(b)=  ^  (ba)^  -  log  b  . 


When  these  expressions  are  used  in  formula  (39)  we  obtain  the  desired 
expression  for  the  likelihood  function: 


A(K,  m  ,  .  .  .  ,  m  ;  a,  b,  u,  a)  =  K 
1  Jx 


1  /s  2  y- 

log  K  -  2  +  —  (bo)  +  log  b 


K 


+  .Z  log  n  (m.-|a,a)  (50) 

1=1  o  1 

In  (50),  of  course,  b  =  1/m  ,  with  m  given  by  (46).  In  the  next  section  we 

e  e 

give  a  numerical  example  of  the  use  of  (50)  and  (46)  to  estimate  |j  and  a  as 
well  as  b. 


A  slight  change  in  viewpoint  to  that  described  at  the  end  of  Section  II 
allows  one  to  invoke  general  properties  of  well-behaved  maximum  likelihood 
estimators.  Specifically,  as  K-ko  our  experimental  data  can  be  considered 
as  resulting  from  larger  and  larger  numbers  of  independent  samples  from 
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a  single  distribution  with  unknown  parameters.  As  such,  the  maximum 
likelihood  estimators  of  these  parameters  (|a,  a,  b)  are  consistent  and  have 
minimum  variance.  This  is,  of  course,  a  much  stronger  statement  than 

A 

our  previous  observation  that  $(b)  $(b)  if  p  and  a  are  known. 

Before  we  proceed  with  numerical  examples,  let  us  see  what  are  some 
further  implications  of  our  assumed  form  for  n(m)  concerning  the  outcome 
of  a  seismicity  experiment.  In  particular,  treating  all  the  parameters  as 
fixed,  let  us  compute  the  expected  number  of  recorded  events,  in  a 

magnitude  interval  m^  <  m  <  Section  II,  the  actual  number, 

^12'  ^  Poisson  variable  (i.e.  ,  (9)  is  valid),  and  the  numbers  of  events 

recorded  in  a  set  of  non-overlapping  magnitude  intervals  are  independent 
random  variables.  From  (30),  on  integrating  by  parts. 


N 


12 


-e 


-  n  (m  -|a,  a)  e 
o  i 


o 


-b(mj  -n) 


(51) 


By  putting  m^  =  m,  m^  =  +oo,  we  obtain  the  expected  number  of  recorded 


events  having  magnitude  at  least  m: 
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—  _  a-bm  ^,rn-|a  ^  ,  a-b|a  + 

N  =  e  erf( - —)  +  e 

m  a 


and  the  expected  total  number  of  events; 


2 

b  a 


,  ^,m-^+ba  ^ 

1  -  erf( - ) 

a 


(52) 


2 

u 

N  =  N  =  ^  2 


(53) 


Given  that  an  experiment  produced  a  total  of  K  recorded  events  (of 
all  magnitudes)  then  the  probability  that  N  ,  the  number  having  magnitude 
at  least  m,  is  equal  to  k  is  given  (just  as  in  (10))  by  the  binomial  distribu¬ 
tion 

Prob  |N^  =  k|N=K|  =  (54) 

where  q  =  1  -  p  and  the  "success  probability,"  p,  is 

p  =  N  /N  .  (55) 

m 

The  conditional  mean  of  N  ,  given  that  N  -  K  is  just  like  (11): 

m  ° 

E  (N  )  =  Kp  (56) 

IX  m 

Equations  (56)  and  (54)  are  used  in  Section  IV  to  indicate  that  the  model 
developed,  including  use  of  the  error  function,  gives  a  reasonable  picture 
of  observed  seismicity  and  the  fluctuations  to  be  expected  in  the  observed 
number  of  events. 
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IV.  EXPERIMENTAL  RESULTS 


We  have  applied  the  maximum  likelihood  method  described  in  the 
preceding  section  to  samples  of  worldwide  data  obtained  from  the  U.S.C. 
and  G.  S.  Our  objective  has  not  been  to  evaluate  the  network  or  to  obtain 
better  estimates  of  worldwide  seismicity.  Rather,  the  objective  has  been 
to  check  the  relevance  to  real  data  of  the  theoretical  machinery  which  has 
been  developed.  It  would  appear  that  the  maximum  likelihood  method  and 
the  statistical  model  used  are  indeed  adequate  and  supply  considerable 
insight  into  the  problem  of  estimating  seismicity. 

The  first  2000  events  of  1968  reported  by  U.S.C.  and  G.  S.  have  been 
used  to  estimate  |a,  a,  a*  and  b*.  The  primed  (base  10)  seismicity  para¬ 
meters  were  estimated  to  be  consistent  with  standard  seismological  usage. 
Constant  terms  were  subtracted  from  the  likelihood  equation  (50)  and  the 
result  has  been  plotted  as  solid  contours  on  Figures  3a  and  3b,  Note  that 
the  contour  intervals  are  not  uniform.  The  function  b*  =  l/(m^  log  10), 
with  m  given  by  (46)  is  also  shown  by  dashed  contour  lines  on  the  figures. 
The  likelihood  function  in  this  example  is  maximized  by  jZi  =  5.  1,  a  =  0.  415, 
and  B*  =  1.725.  The  corresponding  value  of  S  ,  assuming  a  one-year  time 

period,  has  also  been  computed  by  solving  log  K  =  a  -  h[X  H - - —  for  a 

th 

where  K  is  not  2000  but  4500  since  the  2000  event  occurred  on  1  1  June 
1968,  From  this  we  obtained  a*  =  a  log  10  =  11.9.  The  estimated  values  of 
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a*  and  b*  are  in  reasonable  agreement  with  results  obtained  by  Gutenberg 
(3)  using  only  large  events  over  an  extended  time  period.  The  Gutenberg 
values  are  b*  =  1.  59  and  a*  =  12.2  if  it  is  assumed  that  surface  wave  magni¬ 
tudes,  M  ,  are  related  to  m,  by  M  =  1 .  59m  -  3.  97.  (3) 
s  b  s  b 

Figure  4  shows  the  experimental  data  used  in  the  above  experiment. 

The  data  has  been  normalized  to  obtain  the  percent  of  events  observed  above 

magnitude  m  as  a  function  of  m  .  The  theoretical  success  probability 
b  b 

p  =  N  /N,  with  N  given  by  (52)  and  N  =  2000,  has  also  been  shown  for 
m  m 

IJ  =  5.  1 ,  a  =  0.  41  5,  and  b*  =  1 . 725  which  are  the  maximum  likelihood  para¬ 
meter  estimates. 

In  some  applications  of  the  maximum  likelihood  method  of  parameter 
estimation  it  has  been  noted  that  the  estimates  are  also  optimum  in  the 
sense  of  some  reasonable  criterion  for  fitting  data  with  known  parameter¬ 
ized  functions  and  that  the  likelihood  function  is  a  measure  of  the  goodness 
of  fit  as  a  function  of  parameter  values.  The  most  trivial  example  is  the 
estimation  of  the  mean  of  a  Gaussian  distribution  using  N  independent 

samples  from  the  distribution.  In  this  case  the  log  likelihood  function  can 
1  N  2 

be  taken  as  —  .E,  (x.-Li)  where  the  x.  are  the  observations.  Then  the 

N  1=1  1  1 

maximum  likelihood  estimator  is  a  least  mean  square  estimator  and  the 
likelihood  function  is  just  the  mean  square  error. 

We  have  been  unable  to  explicitly  demonstrate  a  simple  cur ve -fitting 
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interpretation  to  the  likelihood  function  for  the  problem  considered  in  this 
paper.  Nevertheless  it  appears  to  be  qualitatively  reasonable  to  accept  the 
assertion  that  the  likelihood  function  is  a  measure  of  the  degree  of  agree¬ 
ment  between  experimental  data  and  the  parameter  values  for  which  the 
likelihood  function  is  evaluated.  For  example,  consider  the  |ji,  a  and  b 
values  which  have  been  indicated  by  the  four  symbols  on  the  -1327  contour 
on  Figure  3b.  A  few  of  the  success  probability  values  implied  by  these 
parameter  values  have  been  indicated  on  Figure  4.  None  of  these  |ji,  a,  and 
b*  values  do  significantly  more  violence  to  the  experimental  data  than  do 
the  maximum  likelihood  estimates.  Thus,  from  the  subjective  viewpoint  of 
fitting  the  experimental  success  probability,  the  small  perturbation  of  the 
likelihood  function  has  had  only  a  small  effect.  Except  for  the  maximum 
likelihood  estimates  no  values  of  success  probability  have  been  shown  below 
m^  =  5.5  since  the  different  values  are  not  clearly  distinguishable  on  the 
figure  at  lower  magnitudes.  Other  values  of  |ji,  a,  and  b*,  which  correspond 
to  much  larger  changes  of  likelihood,  have  also  been  checked  and  found  to 
result  in  large  disagreement  between  the  data  and  theoretical  success 
probability. 

From  the  point  of  view  of  fitting  a  curve  to  the  experimental  data  on 
Figure  4,  it  is  clear  that  significant  changes  in  |jL,  a,  and  b*  from  the  maxi¬ 
mum  likelihood  values  can  have  only  a  limited  effect  if  the  changes  are 
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coordinated  to  have  little  effect  on  the  likelihood  function.  Thus,  moving 
\l  and  a  along  the  ridge  of  A  by  10%  and  making  a  corresponding  change  in 
b*  of  about  30%  has  only  a  limited  effect  on  the  value  of  the  likelihood  func¬ 
tion  and  the  fit  of  the  success  probability  to  experimental  data.  It  should  be 
noted  that  this  discussion  does  not  directly  relate  to  stability  of  the  maxi¬ 
mum  likelihood  estimates.  However,  one  might  conjecture  that  the  shape 
and  steepness  of  the  maximum  of  the  likelihood  function  is  closely  related 
to  stability.  In  this  case,  it  would  seem  that  estimates  of  b*  will  tend  to  be 
less  stable  than  those  for  |a  and  o  if  all  three  are  to  be  estimated.  In 
addition,  the  curvature  of  the  likelihood  function  in  the  region  of  its  maxi¬ 
mum  suggests  that  satisfactory  stability  may  require  long  periods  of 
observation. 

The  stability  of  estimates  is  also  related  to  the  fluctuations  in  experi¬ 
mental  success  probabilities  which  might  be  obtained  for  a  typical  experi¬ 
ment.  Such  fluctuations  have  been  investigated  using  a  short  run  of  U.S.  C. 
and  G.  S.  data  (90  events  from  two  PDE  cards).  The  short  run  of  data  was 
used  to  emphasize  the  fluctuations.  The  likelihood  function  showed  a  very 
broad  minimum,  but  the  parameter  values  b*  =  1. 65,  |a  =  4.  9,  and  a  =  0.  39 
were  near  the  minimizing  set.  The  curve  of  the  success  probability, 
computed  from  equations  (52)  through  (56),  for  these  parameter  values  is 
plotted  as  the  solid  line  on  Figure  5.  The  experimental  data  are  also  shown 
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on  the  figure.  By  using  tables  of  the  cumulative  binomial  distribution  for 
90  trials  and  various  success  probabilities  we  determined  integers  which 
most  nearly  satisfied  the  equations. 


Prob  In  =k/N=K|  =  5% 
;=0  L  m  J 


and 


K 


95 


Prob  |N^  =  k/N=K]-  =  5% 


In  other  words,  and  determine  the  5%  tails  as  used  in  Section  II. 

5  95 

Using  (55)  to  relate  success  probability  to  magnitude  and  normalizing  k^ 
and  k^^  as  a  fraction  of  K,  we  obtain  the  dashed  curves  in  Figure  5.  It 
appears  that  the  model  gives  a  reasonable  picture  both  of  the  average 
seismicity  and  the  fluctuations  to  be  expected  about  the  mean. 
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Fig.  1.  Confidence  intervals  of 
cumulative  histogram  seismicity 
law  given  that  100  events  occurred. 
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Fig.  2.  Confidence  inter¬ 
vals  of  cumulative  histogram 
seismicity  law  given  that  z 

1000  events  occurred. 
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(a) 

Fig,  3.  Lrikelihood  function  contours  and  contours  of 
estimates  of  the  slope  of  the  seismicity  curve. 
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Fig.  3.  Continued. 
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SUCCESS  PROBABILITY 


Fig.  4.  Observed  (2000  U.  S.  C.  and  G.  S.  events)  and 
theoretical  cumulative  histograms. 
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SUCCESS  PROBABILITY 


Fig.  5.  Observed  (90  U.  S.  C.  and  G.  S.  events)  and  theoretical 
cumulative  histograms  and  90%  confidence  intervals. 
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